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Abstract The giant spiny frog (Quasipaa spinosa) is an endangered species with a relatively small distribution limited 
to southern China and Northern Vietnam. This species is becoming increasingly threatened because of over-exploitation 
and habitat degradation. This study provides data on the genetic diversity and population genetic structure of the giant 
spiny frog to facilitate the further development of effective conservation recommendations for this economically 
important but threatened species. We examined 10 species-specific microsatellite loci and Cyt b genes (562 bp) collected 
from 13 wild populations across the entire range of this species. Results of 10 microsatellite loci analysis showed a 
generally high level of genetic diversity. Moreover, the genetic differentiation among all 12 populations was moderate to 
large (overall F;;= 0.1057). A total of 51 haplotypes were identified for Cyt b, which suggests high haplotype nucleotide 
diversities. Phylogeographic and population structure analyses using both DNA markers suggested that the wild giant 
spiny frog can be divided into four distinct major clades, i.e., Northern Vietnam, Western China, Central China, and 
Eastern China. The clades with significant genetic divergence are reproductively isolated, as evidenced by a high 
number of private alleles and strong incidence of failed amplification in microsatellite loci. Our research, coupled with 
other studies, suggests that Q. spinosa might be a species complex within which no detectable morphological variation 
has been revealed. The four phylogenetic clades and some subclades with distinct geographical distribution should be 
regarded as independent management units for conservation purposes. 
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1. Introduction sea level (Zhao, 1998). In Chinese markets, the giant 
spiny frog is especially valued for its medicinal and 


The giant spiny frog (Quasipaa spinosa) is known for the nutritional properties, particularly for the high amount of 


keratinized skin-spines on its chest and the development 
of hypertrophied forearms as a secondary sex organ 
in males during the breeding season (Yu et al., 2010). 
This species inhabits rocky streams in evergreen forests 
and open fields on mountains 500 m to 1,500 m above 
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protein contained in its muscles (Ye et al., 1993; Zhao, 
1998). In addition to harvesting, the habitats of Q. spinosa 
were destroyed by pollution from upland agriculture 
and dam and hydropower projects (Chan et al., 2014). 
Commercial over-exploitation and habitat destruction 
caused the wild population of this species significantly 
decline, both in number and distribution (Zhao, 1998). 

In 2001, the International Union for the Conservation 
of Nature (IUCN) Red List registered this species as 
vulnerable (A2abc) and estimated that over the last three 
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generations (approximately 15 years), the wild population 
had decreased by 30% (Lau et al., 2004; Zhao, 1998). 
Despite the initiation of domestication and artificial 
breeding programs in the 1980s, the culture of giant spiny 
frog remains relatively small in scale because of low 
fertilization and hatching rates, disease, high overwinter 
mortality, and inbreeding depression (Chan et al., 2014; 
Liang et al., 2013). Farming of giant spiny frog mainly 
depends on the natural supply of tadpoles or adults (Chan 
et al., 2014; Liang et al., 2013). 

Recent phylogenetic study suggested that Q. spinosa 
may be a cryptic species complex (Che et al., 2009; Ye 
et al., 2013). However, most frog farms are using frogs 
from different geographical populations to produce 
offspring without a detailed record of genetic and 
performance data (Chan et al., 2014). This practice 
poses a risk of genetic contamination and a threat of 
pathogen transmission. Escapees from aquaculture and 
translocation or even deliberate releases of tadpoles, if 
not sourced from the same location, have the potential 
to interbreed with wild frog, thus possibly modifying 
the genetic structure of wild populations. Interbreeding 
of Q. spinosa from genetically distinct populations can 
potentially reduce the fitness and viability of a species, 
which then becomes increasingly vulnerable to extinction 
over time (Borrell et al., 2011; Miller et al., 2011). These 
issues require knowledge of the genetic structure of the 
species; such knowledge is important for the conservation 
of indigenous genetic diversity (Chan et al., 2014; Miller 
et al., 2011). 

Information on the intra-specific genetic diversity 
and relationships among the wild populations of this 
endangered species remains lacking. Che et al. (2009) 
suggested that the populations currently under Q. spinosa 
likely belong to at least three independent evolutionary 
lineages, but only seven specimens from four localities 
were used as bases from which to infer phylogenetic 
relationships. The sample sizes are too small to reflect 
the genetic diversity and relationships among the 
populations of Q. spinosa. Ye et al. (2013) also supported 
that Q. spinosa is a species complex, rather than a single 
species, on the basis of 12S rRNA and 16S rRNA genes. 
However, the use of mtDNA alone is insufficient to reflect 
the quantification and distribution of genetic variation 
(Dupuis et al., 2012). 

In this study, we analyzed data from microsatellite (10 
microsatellite loci) and mitochondrial (Cyt b gene) loci 
collected from Q. spinosa along the entire range of this 
species to evaluate the phylogeography and population 
genetics of Q. spinosa, as well as to provide genetic 


information for further developing effective conservation 
recommendations to protect this economically important 
but threatened species. 


2. Materials and Methods 


2.1 Sample collection and DNA extraction Two 
hundred and five samples were collected from 13 
localities during 2006 and 2007, in China and Vietnam, 
namely the province (city); Anhui (Huangshan), Fujian 
(Jianyang, Wuyishan), Hunan (Pingjiang), Guangdong 
(Yangshan), Guangxi (Longsheng, Yongfu), Jiangxi 
(Jinggangshan), Yunnan (Pingbian), Zhejiang (Jinhua, 
Lishui) and Northern Vietnam. Specific spatial details 
of the collection sites are given in Fig. 1 and Table 1. 
The collection sites were almostly covered the entire 
geographical range of giant spiny frog (Frost et al., 2006). 
Some samples were released after obtaining their toe 
clips. 

Muscle or toe clips samples were stored in 95% or 
100% ethanol, or frozen at —80°C. The study protocol 
was reviewed and approved by the Committee of Animal 
Research Ethics of Zhejiang Normal University. DNA 
was extracted using the standard phenol-chloroform 
procedure (Sambrook et al., 1989). 


Table 1 Sampling information of giant spiny frog. Geographic 
sampling localities with their coordinates, number of sample sizes 
for sample sizes of Cyt b (N) and number of sample sizes for 


microsatellite (n) are shown. 


Geographic locality (ab.) N/n Coordinates 
Northern Vietnam (VT) 5/0 E 104.33°, N 21.60° 
Pingbian, Yunnan (PB) 5/9 E 103.60° , N 22.81° 
Yongfu, Guangxi (YF) 14/16 E 109.98°, N 24.98° 
Longsheng, Guangxi (LH) 13/27 E 109.98°, N 25.82° 
Yangshan, Guangdong (YS) 4/17 E 112.63°, N 24.48° 
Pingjiang, Hunan (PJ) 16/23 E 113.58°, N 28.72° 
Jinggangshan, Jiangxi (JG) 6/24 E 114.10°, N 26.34° 
Xingan, Jiangxi (XG) 5/10 E 115.60°, N 27.77° 
Wuyishan, Fujian (WY) 8/19 E 118.01°, N 27.27° 
Jianyang, Fujian (JY) 5/13 E 118.60°, N 27.41° 
Lishui, Zhejiang (LS) 11/18 E 119.54°, N 28.27° 
Jinhua, Zhejiang (JH) 15/14 E 119.62°, N 29.11° 
Huangshan, Anhui (HS) 15/15 E 118.57°, N 30.07° 
Total 122/205 
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2.2 Microsatellite genotyping and analysis Among the 
15 developed and characterized species-specific genetic 
markers, 10 proved to be polymorphic and were used in 
this study (Table S1) (Zheng et al., 2009). Microsatellite 
genotyping was performed on all individuals from 13 
localities, but the Northern Vietnam population could not 
be included in the analyses because all individuals and 
loci failed in the amplification. The DNA of Northern 
Vietnam individuals was extracted several times using 
a standard phenol-chloroform procedure and DNeasy 
Tissue Kit (Qiagen), and the electrophoretic band (0.8% 
agarose) was significantly brighter than 25 ng standard 
for DNA quantification. Amplification was performed by 
gradient temperature polymerase chain reaction (PCR) 
(44°C to 68°C), but none of the loci can be successfully 
amplified. 

Amplification was performed in a reaction volume 
of 25 uL with approximately 0.5 uL genomic DNA, 
15.88 uL of ddH,O, 2.5 uL of 10x PCR buffer (Mg” 
free), 1 uL (10 pmol) of each primer, 2 uL MgCl, (25mM), 
2 uL dNTP Mixture (2.5mM), and 0.125 uL of TaKaRa 
Taq (5 U/pL). Following the dye label incorporation 
procedures from Schuelke (Schuelke, 2000), the forward 
primer 5’ was appended with a 19 bp M13 sequence 
(CACGACGTTGTAAAACGAC) and employed at 
0.4 um, whereas the reverse primer used a fluorescent- 
labeled “M13” primer and employed at 0.4 pmol [either 
IRD700 or IRD 800(LI-COR)]. The PCR thermocycler 
profile was as follows: 95°C for 5 min followed by 35 
cycles at 95 °C for 30 s, and annealing temperature for 
30 s (Table S1); 72 °C for 30 s; and a final extension step 
at 72°C for 8 min. The amplified products were analyzed 
on a LI-COR4300 automated DNA sequencer using 6.5% 
denaturing polyacrylamide gels with appropriately labeled 
50 bp to 350 bp internal size standard. The gel images 
were analyzed using the LI-COR SAGA software. 

The number of alleles per locus (Na), heterozygosity 
(H), Allelic frequencies, polymorphic information 
content (PIC), as well as observed (Ho) and expected 
(H) heterozygosities were calculated using Cervus 2.0 
software (Marshall et al., 1998). Unique alleles, defined 
as an allele found in only one of the populations, were 
detected using GDA (Lewis and Zaykin, 2000). The 
potential deviations from Hardy—Weinberg equilibrium 
(HWE) and linkage disequilibrium (LD) were estimated 
using Arlequin version 2.0 and modified Fisher’s exact 
test (Guo and Thompson, 1992; Schneider et al., 2000). 
Results of the multiple comparisons were adjusted using 
Bonferroni correction. 

The population genetic differentiation was quantified 
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using the F,,-estimator of Weir and Cockerham (Weir 
and Cockerham, 1984). A Bayesian clustering method 
was used to organize the dataset and to assign individuals 
to inferred clusters by using the STRUCTURE software 
package (Pritchard et al., 2000). Individuals were placed 
into K clusters using the model-based algorithm, with 
K chosen in advance. A total of 12 independent runs of 
K=1 to K=12 were performed at 1,000,000 Markov Chain 
Monte Carlo (MCMC) repetitions with a 100,000 burn-in 
period assumed with allele frequencies and combinations 
without prior information (Du et al., 2012; Pritchard et 
al., 2000). The log likelihood was used to choose the 
most likely value for K. The posterior probability K 
[PCK|X)] was then calculated. Statistic DK, which is the 
second-order rate of the log probability change of the 
data between successive values of K, was also estimated 
(Evanno et al., 2005). The method successfully detected 
the appropriate number of clusters using multilocus 
genotype data under a number of gene exchange models. 
Given that DK could not be evaluated when K = 1, we 
explored the probability of the K=2 to K=12 on the basis 
of the above results. 

Partitioning genetic diversity through analysis of 
molecular variance (AMOVA) was performed using 
ARLEQUIN version 2.0 on the basis of microsatellite 
data (Schneider et al., 2000). In AMOVA, the total 
variance was divided into variance among groups, among 
populations within groups, among individuals within 
populations, and within individuals. 


2.3 Mitochondrial DNA amplification and analysis 
The Cyt b gene was amplified with primers CB-10933 (5'- 
TAT GTT CTA CCA TGA GGA CAA ATA TC-3’) (Simon 
et al., 1994) and modified CytbC (5'-CTA CTG GTT GTC 
CTC CGA TTC ATG T-3') (Bossuyt and Milinkovitch, 
2000). Amplification was performed in a reaction volume 
of 50 pL that contained approximately 100 ng of template 
DNA, 2.5 mM MgCl,, 2.0 mM 10xPCR buffer, 0.4 uM of 
the forward and reverse primers, 0.2 mM of each dNTP, 
and 1 U rTaq polymerase (TaKaRa, Dalian, China). The 
PCR thermocycler profile was as follows: 94 °C for 2 
min, followed by 35 cycles at 94 °C for 45 s, 56 °C for 45 
s, and 72 °C for | min; and a final extension step at 72 °C 
for 10 min. The PCR products were purified using Wizard 
PCR Preps DNA purification kit (Promega, Beijing, 
China). Each fragment was sequenced with an ABI 
automated DNA sequencer in both directions by using the 
PCR primers. 

All sequences were checked with BioEdit software 
(North Carolina State University) and by visual 
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inspection. The sequences were then imported into 
Clustal X (Thompson et al., 1997) for alignment. Genetic 
diversities for all sampling regions or populations were 
estimated using haplotype (h) and nucleotide diversities 
(x), as implemented in DnaSP version 4.0 (Rozas et al., 
2003). Sequence divergences between haplotypes 
and populations were calculated using Arlequin 
2.0 (Schneider et al., 2000). Pairwise K2P (Kimura 
2-parameter) distances and diagnostic sites for Cyt b 
gene sequences were calculated in MEGA 5.0 to evaluate 
the intrapopulation and interpopulation divergence of 
Q. spinosa (Tamura et al., 2007). The phylogenetic 
relationship among haplotypes was reconstructed using 
maximum likelihood (ML) heuristic search methods 
through by PAUP*4.0 (Swofford, 2002). The robustness 
of these analyses was assessed using non-parametric 
bootstrap replications by PAUP*4.0 with 100 replications 

ML analyses were performed with the general time- 
reversible (GITR+I+G) model of DNA evolution, which 
has the shape parameter of a gamma distribution (G) 
(Posada and Crandall, 1998; Swofford, 2002). The 
GTR+I+G model was identified as the best-fitting model 
under the Akaike information criterion as implemented 
in Modeltest 3.7 (Posada and Crandall, 1998). In the 
Bayesian analysis, we used MrBayes 3.1.2 under the 
GTR+I+G model (Huelsenbeck and Ronquist, 2001). 
Posterior distributions were determined by MCMC 
analysis with one cold chain and three heated chains 
run for 10 million generations with random starting 
trees. Every 100th tree was sampled. Bayesian analyses 
were repeated twice, and the same topology was always 
retrieved. According to the likelihood plots, Inl values 
stabilized with 10,000 generations. Thus, the first 100,000 
generations were discarded as burn-in. Meanwhile, 
Bayesian posterior probabilities were calculated according 
to the remaining set of trees. The rice frog (Fejervarya 
limnocharis) was chosen as the outgroup for all 
phylogenetic reconstructions (Che et al., 2009; Liu et al., 
2005). 

To assess the most probable population configuration 
and geographical subdivision, hierarchical AMOVA was 
performed using Arlequin 2.0 (Schneider et al., 2000). 
The populations were grouped according to different 
geographical distinctions. Genetic differentiation 
among samples from different geographic origins was 
assessed by comparing the average number of pairwise 
differences between populations (PiXY), the average 
number of pairwise differences within populations (PiX 
and PiY), and the corrected average pairwise difference 


(PiXY_(PiX+PiY)/2). Distances and standard error 
were calculated using Arlequin 2.0 with 1000 bootstrap 
replicates (Schneider et al., 2000). 


3. Results 


3.1 Microsatellite DNA analysis 

Genetic variation: All 10 microsatellite loci in the 12 
populations showed high genetic polymorphism (Tables 
S1-S2). However, individuals from Northern Vietnam 
were either non-functional (not amplifying after repeated 
attempts) or poorly amplified (the electrophoretic bands 
of PCR products were dispersed). A total of 339 
different alleles were detected in 205 individuals from 
12 populations. The allele sizes at individual loci 
varied between 12 bp (for locus Psp7) and 92 bp (for 
locus Psp10). Examination of genotyping errors using 
MicroChecker revealed no evidence of large allele 
dropout or stutter-band scoring at any of the 10 loci. 
The mean number of alleles varied from 6.3 to 15.1. The 
mean PIC for the 12 populations ranged between 0.63 
and 0.88. Within all populations, the Ho (0.49-0.72) 
was lower than H, (0.69-0.91) (Table S2). The 12 
populations showed very high gene diversity (mean 
He=0.82). Among the populations, the gene diversities 
of XG (He=0.90) and WY (He=0.91) were higher than 
those of other populations. Deviation from HWE and 
linkage disequilibrium (LD) was also evaluated in all 12 
populations. Out of 120 loci (12 populations multiplied 
by 10 loci), 65 showed deviation from HWE (P < 0.05; 
Table S2). All 10 microsatellites were not in LD (P > 
0.05) for each pair of loci across all samples. 


Genetic structure and gene flow: After Bonferroni 
correction, some multi-loci F';; values were small, but 
all values were significant at P < 0.05 (Table 2). All 
of the F,,; values, except for five lower than 0.05, fell 
between 0.05 and 0.25, indicating that a moderate or large 
differentiation existed in the populations (Wright, 1978). 
The average of F,r of all loci was 0.1057. The PJ and HS 
populations had the largest amount of differentiation at 
0.2095. 

The clustering method revealed that AK max was 3.02 
and peaked at K = 3. The highest posterior probability 
(In Pr (X|K) = —4443.6) of the 20 runs for K = 3 was 
presented along with color-coded probabilities of 
individuals belonging to a cluster (Figure 1). The highly 
non-random association of colors showed three distinct 
genetic cluster groups among the sampling localities. 
The cluster groups are the Western China, Central 
China, and Eastern China clusters. The cluster names 


No. 2 Danna YU et al. 
refer to the sample range. The Western China cluster 
comprised alleles from three populations sampled from 
Pingbian (PB) of YUNNAN Province, Yangshan (YS) of 
Guangdong Province, and Longsheng (LH) of Guangxi 
Zhuang Autonomous Region. The Central China cluster 
comprised alleles from populations at Yongfu (YF) of 
Guangxi Zhuang Autonomous Region, Pingjiang (PJ) of 
Hunan Province, Jinggangshan (JG) and Xingan (XG) 
of Jiangxi Province, and Wuyishan (WY) of Fujian 
Province. The eastern China cluster comprised alleles 
from populations at Jianyang (JY) of Fujian Province, 
Lishui (LS) and Jinhua (JH) of Zhejiang Province, and 
Huangshan (HS) of Anhui Province. 

A total of 118 unique alleles were detected among the 


Genetic Diversity and Population Structure of Giant Spiny Frog 79 


three clusters. The Central China cluster possessed the 
highest number of unique alleles (57.63%), followed by 
Western China (23.73%) and Eastern China (18.64%) 
(Table S3). The cluster groups were also analyzed 
using hierarchical AMOVA method, which revealed the 
following statistical significant (P < 0.001) substructure: 
66.01% of the total genetic variations was observed 
among the individuals, 1.95% among the groups, 9.46% 
among populations within the groups, and 22.58% among 
individuals within groups (Table 3). 


3.2 Mitochondrial DNA analysis 

Sequence variation and genetic distance: A total of 218 
variable (38.79%) sites and 211 parsimony-informative 
(37.54%) sites were found in the 562-bp Cyt b fragments. 


Table 2 Pair-wise F,, between the twelve populations estimated from the ten microsatellite loci. 


HS JH LS JY PJ XG WY JG YS LH YF 
HS 
JH 0.2052 
LS 0.1942 0.1299 
JY 0.1673 0.1268 0.0784 
PJ 0.2095 0.1758 0.1310 0.1306 
XG 0.1571 0.1295 0.0808 0.0898 0.0891 
WY 0.1187 0.1213 0.0909 0.0515 0.0942 0.0379 
JG 0.1482 0.1414 0.1118 0.1001 0.1254 0.0937 0.0636 
YS 0.1461 0.1153 0.0669 0.0726 0.1005 0.0615 0.0494 0.0603 
LH 0.1623 0.1328 0.1167 0.1082 0.1274 0.0937 0.0681 0.0848 0.0847 
YF 0.1616 0.1488 0.1313 0.1066 0.1396 0.0902 0.0695 0.1002 0.0882 0.0287 
PB 0.1368 0.1123 0.0904 0.0881 0.1097 0.0686 0.0359 0.0604 0.0529 0.0431 0.0697 


Table 3 Nested Analysis of Molecular Variance (AMOVA) indicating partitioning of genetic diversity. 


Source of Variation d.f. Sum of Squares Variance Components Percentage of Variation 
Among groups 2 66.169 0.090Va 1.95" 
Among populations within groups 9 176.139 0.436 Vb 9.46" 
Among individuals within populations 193 987.793 1.040 Vc 22.58" 
Within individuals 205 623.000 3.039 Vd 66.01" 

Total 409 1853.100 4.604 


“All indices are highly significant; the probability (P) of obtaining a more extreme component estimate by chance alone = <0.001 (using 1000 


randomizations). 
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Figure 1 Color scheme: Bayesian cluster analysis of the microsatellite variation among the 12 giant spiny frog populations. Each vertical 
line represents an individual. The color composition displays the probability of belonging to each of the three clusters as defined by 
STRUCTURE. The three colors (red, blue, and green) represent the three genetic clusters. Map: Map of sampling localities for this study. 
The 13 sampled localities have colors corresponding to the genetic cluster wherein the majority of their respective individuals were assigned. 
Table 1 shows the abbreviation and coordinates of localities. The dotted area shows the distribution of giant spiny frog. 


These sites defined a total of 51 haplotypes (Figure 2). 
The haplotypic sequences were deposited in GenBank 
with accession Nos. EU156068—156118. Each population 
has two to seven haplotypes. Each population also 
contained unique haplotypes. 

Haplotype and nucleotide diversities were respectively 
h=0.968+0.006 and z=0.1152+0.0068 for all samples 
combined. The nucleotide diversity had an overall trend 
of decline from west to east (Table 4). Estimates of 
sequence divergence (K2P) among populations ranged 
from 0.419% (Jianyang, Fujian Province vs. Longsheng, 
Guangdong Province) to 24.37% (Northern Vietnam vs. 
Wuyishan, Fujian Province) (Table S4). 


Phylogenetic relationship and population structure: 
The ML and Bayesian inference produced highly similar 


topologies, which clustered the haplotypes into four major 
clades, i.e., Northern Vietnam, Western China, Central 
China, and Eastern China clades (Figure 2). The clade 
names refer to their geographical pattern. The Northern 
Vietnam clade has three unique haplotypes, which formed 
the most basal branch in all phylogenetic reconstructions. 
Haplotypes H39 and H41 were shared in the Central and 
Eastern China clades despite having a closer relationship 
with other haplotypes of the Eastern China clade. Genetic 
distances among four major clades (range: 7.121—20.027) 
were almost one order of magnitude higher than those 
within each group (range: 0.925-8.671; Table 4). To 
provide a statistical evaluation of the geographical 
patterns suggested by the phylogenies, all possible sets 
of AMOVA comparisons were calculated. In the schemes 
dividing the giant spiny frog populations into two, three, 
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Figure 2 Phylogenetic relationships of haplotypes, and geographic distribution and frequency of each haplotype in the giant spiny frog. 
Phylogenetic relationships of haplotypes based on the Cyt b gene sequences were derived from Maximum-likelihood (GTR+I+G model) 
algorithm and Bayesian inference. The numbers above and below the branches are bootstraps and posterior probabilities values, which only 
shows values higher than 50%. Fejervarya limnocharis was the outgroup species for rooting the tree. 


four, and five groups, the grouping scheme of four groups 
had a high Per value (Øer = 0.0715, P < 0.01; Table S5). 
Genetic distances among the four major clades of giant 
spiny frogs ranged from 7.12% (Eastern China clade vs. 
Central China clade) to 20.03% (Northern Vietnam clade 
vs. Central China clade) (Table 5). 


Diagnostic sites: Figure S1 shows the 30 diagnostic 
sites that clearly distinguished the four major clades. The 
special base composition at sites 41, 164, 458, and 463 
differentiated the alleles H34—51 from the Eastern China 
clade with all alleles from other three major clades, and 
the “C” at site 98 only occurred in the haplotypes from 
the Central China clade (Figure S1). Diagnostic sites were 
also found between certain subclades, e.g., subclade W1 
vs. subclade W2 and subclade E1 vs. subclade E2 (Figure 
S1). 


4. Discussion 


4.1 Genetic diversity within Q. spinosa This study is 
the first to analyze the genetic diversity and structure of 
wild giant spiny frog populations using a multiplex set of 
microsatellite markers. Our result showed that the PIC of 
the majority of the markers was high (92%, range: 0.23- 
0.94), which indicated that the 10 microsatellite markers 
were highly polymorphic and can thus be used for genetic 
diversity and population structure analyses of giant spiny 
frog (Botstein et al., 1980). The average allele number 
of the 10 loci in giant spiny frog was 33.9/locus, which 
is significantly higher than the average allele number in 
some Anura species (Andersen et al., 2004; Burns et al., 
2004). The allele number is usually positively associated 
with sample size (Goldstein et al., 2005). However, the 
high polymorphism of the 10 microsatellites enabled 
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Table 4 The haplotype diversity (h) and nucleotide diversities (z) 
of Cyt b gene in each population of giant spiny frog. 


Location h m 

Northern Vietnam (VT) 0.800+0.164 0.00925+0.00248 
Pingbian (PB) 0.400+0.237 0.09181+0.05447 
Yongfu (YF) 0.758+0.084 0.02951+0.01893 
Longsheng (LH) 0.84840.074 0.09005+0.01829 
Yangshan (YS) 1.000+0.126 0.09401+0.03056 
Pingjiang (PJ) 0.825+0.071 0.00340+0.00058 
Jinggangshan (JG) 0.905+0.103 0.09929+0.02933 
Xingan (XG) 1.000+0.126 0.04057+0.01753 
Wuyishan (WY) 0.750+0.063 0.02440+0.01413 
Jianyang (JY) 0.800+0.164 0.00321+0.00069 
Lishui (LS) 0.727+0.144 0.00285+0.00100 
Jinhua (JH) 0.25740.142 0.00047+0.00027 
Huangshan (HS) 0.752+0.056 0.00275+0.00026 
Total 0.968+0.006 0.11520+0.00680 


the researchers to use powerful statistical tools to 
analyze genetic diversity and population structure in 
12 populations. HWE is an important prerequisite for 
population genetic analysis, but the observed allele 
frequency in more than 50% of the loci in each population 
was found significantly deviate from HWE in some 
studies (Musammilu et al., 2014; Zhou et al., 2004). The 
significant deviation is generally caused by insufficient 
sample size, null alleles, loss of random mating, and the 
Wahlund effect (Rousset and Raymond, 1995; Wahlund, 
1928; You et al., 2008). AMOVA of microsatellites also 
revealed that variation among populations was 1.59%. 
Therefore, the Wahlund effect can be excluded because 
of low genetic differentiation among populations. We 


detected null alleles in four loci (Psp6, Psp9, Psp10, and 
Psp14). However, low frequencies (0.08—0.15) of null 
alleles cannot result in genetic disequilibrium. Non- 
amplifying samples in repeated trials were not present 
in any of the primer pairs in the 12 populations of Q. 
spinosa. The deviation of the loci in this study may 
be attributed to the relatively small sample size and 
heterozygote deficiencies, which are caused by selection, 
population mixing, or nonrandom mating (Rousset and 
Raymond, 1995). 

As shown by the significant departure from HWE at 
some microsatellite loci, the observed heterozygosities 
were relatively low in all populations (Table S2). A 
heterozygote deficit is likely to be apparent in Q. spinosa 
from the 12 populations, and heterozygosity might be 
more vulnerable to selection pressures and/or habitat 
changes (Musammilu et al., 2014). No historical data 
were available on the genetic diversity and heterozygosity 
of Q. spinosa but departure of the microsatellite data from 
HWE expectation in all populations may be attributed 
to commercial over-exploitation and habitat destruction. 
Although Q. spinosa is listed as vulnerable by the IUCN 
and the China Red List (Lau et al., 2004; Zhao, 1998), 
many populations accessible to hunters may already have 
been extirpated because such populations may be hunted 
multiple times per year by multiple hunters (Chan et al., 
2014). 

Despite using a similar method of analysis, the mean 
Ho levels found in this study (0.49-0.72) are higher than 
those found by Andersen et al. (2004) for the Danish 
European tree frog populations (0.35-0.53) and by 
Newman and Squire (Newman and Squire, 2001) for the 
common and widespread Rana sylvatica (0.44—0.50). 
However, this results of this study are similar to those on 
the formerly widespread Litoria aurea (0.43-0.82) by 
Burns et al. (2004). Higher genetic variation in wild giant 
spiny frog populations can positively affect the protection 
of giant spiny frog. 


Table 5 Pairwise genetic distances for Cyt b sequences of the giant spiny frog. Below diagonal: Net genetic distances between the four major 


phylogenetic clades [PiXY—(PiX+PiY)/2]. Above diagonal: their standard errors of distances between the four major phylogenetic clades. 


Values within clades and their standard errors in parenthesis are shown in the diagonal (PiX). 


Northern Vietnam 


Western China 


Central China Eastern China 


Northern Vietnam 0.925(0.248) 2.441 2.267 2.266 
Western China 17.245 8.671(1.384) 1.134 1.421 
Central China 20.027 13.523 5.186(1.029) 0.651 
Eastern China 18.451 14.176 7.121 1.566(0.078) 


No. 2 Danna YU et al. 


4.2 Phylogeography and population structure 
Amphibians, especially frogs, have species recognition 
and mate choice systems that rely on non-morphological 
characteristics (e.g. advertisement calls) and have 
the tendency to exhibit conservative morphological 
evolution. Thus, frogs are considered to have a large 
amount of cryptic genetic diversity (Bickford et al., 
2007; Funk et al., 2012; Stuart et al., 2006). In the 
Vietnam population, the microsatellite loci were either 
non-amplifying (Vietnam) or have a high number of 
private alleles in clades (Table S3). Differentiation at 
microsatellite loci among clades can indirectly support 
reproductive isolation. Q. spinosa exhibits high site 
fidelity, with some individuals found in similar locations 
in successive years (Chan et al., 2014; Palo et al., 2004). 
The microsatellite dataset on the 12 giant spiny frog 
populations (except Vietnam) exhibited a moderate or 
large genetic differentiation with an overall F, of 0.1057, 
which was expected for the frog with low dispersal rates. 
Q. spinosa is thus expected to exhibit a higher degree of 
population subdivision than any other major animal taxa 
because of poor dispersal capabilities (Chan et al., 2014; 
Palo et al., 2004). STRUCTURE suggested that the frog 
populations can be split into three distinct genetic clusters, 
i.e., Western China, Central China, and Eastern China, 
respectively (Figure 2). Although AMOVA only explained 
1.95% of the observed genetic variation among groups, 
it was statistically significant (P<0.001). The genetic 
variation confirmed the existence of differentiation 
among the three groups. The phylogenetic trees using 
Cyt b gene further corroborated that the four clades are 
evolutionarily distinct (Figure 2). Relationships among 
the samples from different geographical locations and 
putative phylogeographic groups using microsatellite and 
mitochondrial markers were congruent with the results of 
Ye et al. (2013). 

The clustering result matched the geographic 
distribution and geological topography of the populations. 
Giant spiny frogs distributed in the eastern areas are 
isolated from the frogs in central regions by strong 
biogeographic barriers, such as Poyang Lake, Gan River, 
and Wuyi Mountains. The phylogenetic reconstructions 
based on Cyt b gene also generated a topology that 
showed a stepwise-ladder pattern (Figure 2). Out of the 
four clades, the Northern Vietnam clade was the most 
basal branch, followed by the divergence of Western 
China clade and the Central and Eastern China clades. 
The west-to-east dispersal can also be corroborated by the 
gradual decrease in genetic diversity from west to east. 
Populations from the original location generally usually 
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lose their genetic diversities during dispersal or expansion 
(e.g., Yang et al., 2004). 

The large number of variable sites (38.79%) and 
parsimony-informative sites (37.54%) in Cyt b gene also 
proved that Q. spinosa may be a species complex. Genetic 
distances among the four major clades of giant spiny 
frogs are higher than or comparable with those reported 
in previous studies (Barber, 1999; Bradley Shaffer et al., 
2004; Elmer et al., 2007; Yu et al., 2015), as shown by 
interspecific differentiation in most other frogs. Barber 
(1999) found mean pairwise distances ranging from 
12% to 20% for mitochondrial Cyt b sequences among 
different Hyla species (i.e., H. arenicolor and H. eximia), 
which are comparable with the divergences among the 
four major clades in this study. Some diagnostic sites 
and unique alleles were also found to distinguish each 
subclade pair reliably (Figure S1; Table S3). 

However, STRUCTURE showed that the YF and YS 
populations were respectively clustered to the Central and 
Western China clades. This result differed from those of 
mitochondrial analysis (Figure 1). Admixture structure 
also existed in some populations, i.e., YF population 
mixed with the cluster of Eastern China despite 
being near the cluster of Western China according to 
geographical distance. The XG and WY populations also 
mixed among the three clusters (Figure 1). Discordance 
between mtDNA and nDNA genes can be attributed 
to introgression between distinct groups, differences 
between male and female dispersal rates, or incomplete 
lineage sorting of ancestral polymorphism (Qu et al., 
2012; Yang and Kenagy, 2009). 

Our data rejected the possibility that introgression 
caused the pattern. Introgression always occurred in 
contact-zone individuals (Yang and Kenagy, 2009). 
However, admixture was found among far-distance 
populations, i.e., XG and WY populations showed a 
mixing among three clusters. Another three populations 
from Central China only mixed with the cluster of Eastern 
China. The population of Western China (e.g., LH) mixed 
with the cluster of Eastern China. As an alternative, we 
considered that the retention of ancestral polymorphism 
at the nuclear microsatellite loci likely caused the mixing. 
The uniparental inheritance of mtDNA resulted in an 
effective population size four times smaller than that of 
bi-parentally inherited nuclear DNA. Thus, nuclear DNA 
is more likely to retain ancestral polymorphism (Yang and 
Kenagy, 2009). However, further systematic sampling 
from more areas and more molecular markers should 
be used in future investigations to confirm the present 
hypothesis. 
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4.3 Implications for conservation Our research coupled 
with other studies suggests that Q. spinosa might be 
a species complex (Che et al., 2009; Ye et al., 2013). 
The failure to diagnose biological diversity can hamper 
conservation efforts and basic scientific inquiry (Mayden 
and Wood, 1995). Conservation decisions are often 
made on the assumption that named taxonomic units 
represent evolutionary lineages. Phylogeographic analysis 
provides valuable information on how intraspecific 
genetic variation is partitioned. Such analysis also 
identifies evolutionary significant units (ESUs) and 
management units (MUs) for species. According to the 
criterion proposed by Moritz (1994), populations that are 
differentiated for mtDNA haplotypes or at nuclear loci are 
designated as MUs, whereas those that are reciprocally 
monophyletic for mtDNA haplotypes and show 
significant differentiation at nuclear loci are given ESU 
status. The molecular data presented in this work indicate 
significant geographical structuring of genetic variation 
across the range of Q. spinosa, and the topologies of all 
phylogenetic trees derived from all sequences (Figure 
2) showed that reciprocal monophyly criterion for ESU 
designation was met by the Northern Vietnam, Western 
China, Central China, and Eastern China clades. 

A precise and correct delimitation of species is 
essential because species are basic units for analysis in 
biogeography, ecology, macroevolution, biodiversity 
assessment, and conservation and management (Francesca 
et al., 2006; Sites Jr and Marshall, 2004). Over- or under- 
resolving species boundaries can result in erroneous 
decisions for management and conservation (Francesca et 
al., 2006; Sites Jr and Marshall, 2004). This suggests that 
a re-evaluation of Q. spinosa boundaries and conservation 
status is needed. Our results show that Q. spinosa can 
be divided into multiple ESUs and all of them should be 
considered as evolutionary independent units, and hence 
should be independently managed wherever possible to 
avoid the break-up of coadapted gene complexes and 
other forms of outbreeding depression (Tallmon et al., 
2004; Ficetola and De Bernardi, 2005; Ficetola et al., 
2007). At present, most frog farms are using Q. spinosa 
from different geographical populations, which will lead 
to arisk of genetic contamination and a threat of pathogen 
transmission. So, it is necessary to make periodic genetic 
evaluations of the source population to ensure that they 
are genetically similar (Borrell et al., 2011). 

To protect the relatively abundant genetic variations 
held in different clades, the four phylogenetic clades and 
some subclades with distinct geographic distribution 
may be regarded as independent management units for 


conservation purposes. 
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Table S1 Characteristics of the ten polymorphic microsatellite loci isolated from the giant spiny frog. 


: i — Size range Ta 

Locus Repeat Motif Primer Sequences (5'— 3") (bp) (°C) Source 

Psp1 (GT)is F: 5'-CATTGGCACTGCTGTTATCC-3' 213-243 59 Zheng et al. (2009) 
R: 5'-GACCTGTGACAGCTTCTTTGTG-3’ 

Psp2 (GT)7T(TG)e F: 5'-AACAGTGAAAGAACCGAAAC-3' 142-178 51.5 Zheng et al. (2009) 
R: 5'-CCCACAATGGAATGGACACG-3’ 

Psp6 (GT)i2...(GT)A F: 5-CACGCAAAGTGAAACGCT-3' 291-351 65 Zheng et al. (2009) 

T(GT)e R: 5'-CTCCCCAACGACACACG-3' 

Psp7 (GT)i3 F: 5’-ACTCTACAGCAAGTAAAAGC-3’ 161 - 173 52.5 Zheng et al. (2009) 
R: 5'-CCCAAAATCAGAAAATAAT-3' 

Psp8 (CA)? F: 5'- TGCTTGGTAGTTTGCGATT-3' 314-358 62 Zheng et al. (2009) 
R: 5'- CGTGACCGGAGTGATGTC-3' 

Psp9 (AC) F: 5-GGATCAGCAAGCAACATTAT-3' 220-236 49 Zheng et al. (2009) 
R: 5'-CAGTCCCAGTTTCTTCTCAC-3' 

Psp10  (CG)(CA)4s F: 5- AAAACACAAACAAAAGAA-3' 272-364 55 Zheng et al. (2009) 
R: 5'- AACAGGGTAAATATGGAC-3' 

Psp11  (GT)s F: 5-GGACAGGGTGAAGGCAGTAT-3' 224-274 55,9 Zheng et al. (2009) 
R: 5'-CCTGTGAGGCAATATGAAAA-3' 

Psp12  (GT)x»7 F: 5-TTAAAATAGCCAACCAG-3' 234-286 50 Zheng et al. (2009) 
R: 5'-CATCAATTACCACATGC-3’ 

Psp14 (TG)11 F: 5'-ATGGCTGGTGGAAAAAGACT-3' 228-300 59 Zheng et al. (2009) 
R: 5-TAGGAGGGGCAACGGAG-3' 


Table S2 Allele number (Na), observed (Ho) and expected (He) heterozygosity and 


polymorphic information content (PIC) in the twelve populations of giant spiny frog. 


Pop Parameters san 
Psp8 Psp6 Psp7 Psp2 Psp9 Psp14 Psp10 Psp1 Psp12 Psp11 Mean 
HS Na 2 6 10 11 8 7 6 6 5 2 6.3 
Ho 0.00 0.29* 0.60* 0.73 0.87 0.80 0.93 0.47 0.27 1.00* 0.6 
He 0.27 0.69 0.87 0.88 0.83 0.75 0.79 0.77 0.54 0.52 0.69 
PIC 0.23 062 0.82 0.84 0.78 0.70 0.73 0.71 O49 0.38 0.63 
JH Na 6 10 4 4 5 9 6 9 7 3 6.3 
Ho 0.64* 0.57 0.29* 0.14* 0.43* 0.71 0.50 064 0.50 0.43 0.49 
He 0.72 0.84 0.55 0.56 0.62 0.85 0.81 089 0.78 0.46 0.71 
PIC 0.65 0.79 049 046 0.54 0.80 0.74 0.84 0.72 0.40 0.64 
LS Na 9 10 9 11 3 14 6 14 10 12 9.8 
Ho 0.78 0.28* 0.83 0.83 0.28* 0.78 0.44 0.67 0.61* 0.61 0.61 
He 0.76 0.88 0.80 085 0.56 0.92 0.66 088 0.86 0.79 0.8 
PIC 0.71 084 0.75 O81 045 0.88 0.61 0.85 0.81 0.75 0.75 
JY Na 10 10 8 8 5 12 9 12 9 14 9.7 
Ho 0.85 0.38* 0.77 0.31 0.77 0.54* 0.69 0.62* 0.69 £0.85 0.65 
He 0.89 0.87 0.73 0.53 0.74 0.90 0.85 0.90 088 0.94 0.82 
PIC 0.84 082 0.68 049 0.66 0.85 0.80 086 082 0.90 0.77 
PJ Na 7 17 9 6 3 19 18 14 10 15 11.8 
Ho 0.48 0.74 0.61 O65 0.13* 0.87 0.57* 0.52* 0.57 0.43* 0.56 
He 0.53 0.94 0.70 0.70 0.57 0.94 0.94 092 085 0.84 0.79 
PIC 0.50 091 0.67 0.63 047 0.92 0.92 0.90 082 0.81 0.76 
XG Na 10 11 11 10 5 10 15 9 11 9 10.1 
Ho 0.80* 0.40* 1.00 0.50 0.60* 0.40* 0.70* 0.40* 0.70 0.50* 0.60 
He 0.92 094 0.89 0.89 0.76 0.93 0.97 0.91 0.91 0.89 0.90 
PIC 0.86 089 084 083 068 O87 0.92 085 085 £0.83 0.84 
JG Na 10 12 11 11 6 18 17 16 26 15 14.2 
Ho 0.58* 0.79 0.38* 0.54* 0.33 0.67 0.54* 0.63* 0.75* 0.83 0.60 
He 0.78 089 0.65 085 0.40 0.93 0.92 0.93 0.97 0.93 0.83 
PIC 0.75 0.86 0.61 0.81 O37 0.91 0.90 0.90 0.94 0.90 0.80 
WY Na 15 14 15 11 7 18 18 17 17 19 15.1 
Ho 0.74 0.32* 0.63* 0.63 0.58 0.53* 0.53* 0.68 0.56* 0.58* 0.58 
He 0.89 088 0.90 0.91 0.83 0.95 0.95 0.94 0.92 0.94 0.91 
PIC 0.86 085 0.87 088 0.79 0.92 0.92 091 088 0.91 0.88 
LH Na 9 10 12 15 7 12 17 18 15 12 12.7 
Ho 0.93* 0.58* 0.85 0.56* 0.48 0.37* 0.56* 0.68% 0.41* 0.70 0.61 
He 0.75 0.90 084 092 056 0.87 0.90 093 091 0.83 £0.84 
PIC 0.71 0.87 O81 0.90 050 0.84 0.88 0.90 088 0.79 0.81 
YF Na 7 10 8 11 6 7 10 9 14 9 9.1 
Ho 1.00 0.69 0.81 0.75 0.38 0.56 0.63 0.56* 0.56* 0.63* 0.66 
He 0.77 0.90 0.66 0.92 0.72 0.83 0.91 088 0.92 0.89 0.84 
PIC 0.71 086 0.62 088 066 0.77 0.87 O83 O88 £0.85 0.79 
YS Na 14 12 10 9 4 13 16 14 19 10 12.1 
Ho 1.00 0.59 0.76 0.47* 0.59 0.76 0.76* 0.76* 0.76* 0.76 0.72 
He 0.88 080 0.84 0.79 0.64 0.90 0.94 0.93 0.96 0.90 0.86 
PIC 0.84 0.76 0.79 0.74 0.57 0.86 0.90 090 093 0.86 0.82 
PB Na 11 8 9 9 7 16 10 10 11 8 9.9 
Ho 0.71* 0.29* 0.64* 0.71 0.71 0.71* 0.79 0.57* 0.86 0.79 0.68 
He 0.92 0.87 0.88 0.83 065 0.96 0.86 090 091 0.85 0.86 
PIC 0.87 0.82 0.83 0.78 0.60 0.92 0.81 0.85 0.87 0.80 0.82 


An asterisk (*) indicates the locus which showed deviation from Hardy-Weinberg equilibrium. 


Table S3 Numbers of unique alleles at the ten microsatellite loci in the three giant spiny frog groups. 


Locus 
Clusters Total 
Psp1 Psp2 Psp6 Psp7 Psp8 Psp9 Psp10 Psp11 Pspl2  Psp14 
Eastern China 1 2 7 3 0 1 3 1 2 2 22 
Central China 4 2 6 4 6 4 9 11 12 7 65 


Western China 2 4 1 0 4 0 1 3 4 2 21 


Table S4 


Mean genetic distance between population (below diagonal) and within group divergences (on the diagonal) based on K2P model 
from Cyt b data in the giant spiny frog. 


Population HS JH LS JY WY XG JG PJ YS LH YF PB VT 
HS 0.00275 
JH 0.02784 0.00047 
LS 0.02521 0.01365 0.00285 
JY 0.02568 0.01314 0.00419 0.00321 
WY 0.07873 0.07415 0.07462 0.07460 0.02440 
XG 0.07158 0.06705 0.06566 0.06633 0.03404 0.04057 
JG 0.10545 0.11199 0.11280 0.11418 0.09687 0.09707 0.0839 
PJ 0.07075 0.07745 0.07919 0.07934 0.05376 0.05173 0.07914 0.00343 
YS 0.08612 0.08153 0.07482 0.07712 0.11461 0.10964 0.12436 0.11826 0.09151 
LH 0.15887 0.14241 0.16236 0.16475 0.14945 0.14964 0.15391 0.14137 0.16424 0.08005 
YF 0.19388 0.19623 0.19530 0.19861 0.18755 0.18625 0.18259 0.18436 0.17875 0.07165 0.02951 
PB 0.19309 0.19645 0.19480 0.19763 0.18588 0.18583 0.17950 0.18543 0.17417 0.17977 0.17644 0.07181 
VT 0.21151 0.21679 0.21625 0.21691 0.24371 0.24007 0.22821 0.22698 0.21245 0.21969 0.21030 0.18424 0.00925 


Table S5 Hierarchical analysis of AMOVA of the giant spiny frog from Cyt b data. 


Among pops within Within pops ®sr 
Groups 


groups ®sc 

2 Groups [HS, JH, LS, JY, WY, XG, JG, Saas rea 
PJI[YS, LH, YF, PB, VT] 
3 Groups [HS, JH, LS, JY, WY, XG, JG, 
PIES LH, YE PEINT 0.0891 0.3405 
4 Groups [HS, JH, LS, JYI[WY, XG, JG, 
PJ, YS] LH, YF, PBJ[VT] 0,0870 Paal 
5 Groups [HS, JH, LS, JYI[WY, XG, JG, 
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